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, .  ABSTRACT 

V 

We  consider  the  steady  thermocapitlary  motion  in  a  square  cavity  with  a  top  free  surface  in  the 

absence  of  gravitational  forces.  The  cavity  is  heated  from  the  side  with  the  vertical  boundaries 

isothermal  while  the  horizontal  boundaries  are  adiabatic.  The  relative  change  in  the  surface  tension  is 

very  small,  i.e.  an  appropriate  capillary  number  tends  to  zero,  so  that  the  free  surface  is  assumed  to 

remain  flat  at  leading  order.  A  finite  difference  method  is  employed  to  compute  the  flow  field. 

Numerically  accurate  solutions  are  obtained  for  a  range  of  Prandtl  numbers  and  for  Reynolds 
f QOQO 

numbers  as  high  as  5  x  Surface  deflections  are  computed  as  a  domain  perturbation  for  small 
capillary  number.  In  addition,  asymptotic  methods  are  used  to  infer  the  boundary  layer  structure  in 
the  cavity,  in  the  limit  of  large  values  of  the  Reynolds  number  Re  and  the  Marangoni  number  Ma^For 
a  fixed  Prandtl  number  Pr,  it  is  shown  that  the  Nusselt  number,  liquid  circulation  and  maximum 
vorticity  are  asymptotic  to  Re173,  Re'1/3  and  Re273,  respectively.  These  results  are  in  agreement  with 
the  computed  solutions.  The  leading  order  solution  for  the  free  surface  deformation  is  sensitive  to  the 
value  of  Pr.  With  Pr  >  t ,  the  depression  near  the  hot  corner  may  exceed  the  elevation  near  the  cold 
comer.  While  a  secondary  elevation  may  be  induced  near  the  hot  comer  when  Pr  <  1 . 
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1.  Introduction 

Convective  motions  driven  by  a  temperature  gradient  along  the  interface  between  two  immiscible 
fluids,  due  to  the  variation  of  surface  tension  with  temperature,  are  of  considerable  interest  and  play 
an  important  role  in  small  scale  and/or  low  gravity  hydrodynamics  (Schwabe,  1981  and  Ostrach, 
1982).  Because  these  thermocapillary  flows  occur  in  crystal  growth  melts  and  dominate  the 
convective  flows  in  the  microgravity  environment  of  space,  there  have  been  a  number  of  recent 
studies  of  simplified  two-dimensional  models  with  negligible  gravitational  effects. 

Sen  &  Davis  (1982)  consider  steady  thermocapillary  convection  in  a  differentially  heated 
rectangular  slot  with  a  top  free-surface.  They  present  results  valid  for  vanishingly  small  aspect  ratio  A 
(height/width)  with  the  Reynolds  number  Re  ~  0(A),  the  Marangoni  number  Ma  —  O(A)  and  the 
capillary  number  Ca  —  0(A4).  The  same  problem  is  also  studied  by  Strani,  Piva  &  Graziani  (1983)  for 
A  — *  0  with  Ca  —  0(A4)  but  with  milder  restrictions  on  Re  and  Ma.  These  studies  are  instructive  in 
giving  the  flow  field  and  surface  deflection  in  a  conduction-dominated  regime,  but  give  little 
information  about  strongly  convective  flows. 

Strani  et.  al.  (1983)  also  compute  the  thermocapillary  motions  in  a  rectangular  cavity  with  A  =  0.2,  1 
and  5  and  with  the  Prandtl  number  Pr  =  1  (i.e.  Ma  =  Re).  They  allow  for  free-surface  deformation 
and  conclude,  as  expected,  that  for  Ca  =  0.1  surface  deformation  have  a  negligible  influence  on  the 
flow  field.  Because  they  use  a  coarse  mesh  for  their  finite-difference  computations,  their  results  are 
accurate  only  for  low  Re  and  Ma.  Since  vigorous  motion  occurs  in  float-zones  at  high  values  of  Ma 
and  Re  (Schwabe  1981,  Ostrach  1982),  there  is  a  need  for  further  study  of  this  model  problem. 
Axisymmetric  modeling  of  a  half-floating  zone  configuration  has  been  considered  by  Fu  &  Ostrach 
(1983)  where  they  compute  the  flow  field  at  different  values  of  Re,  Ma  and  Pr.  Many  of  the  features  of 
their  solutions  are  common  to  oursi 

The  stability  of  such  flows  is  of  great  interest.  As  the  experiments  of  Preisser,  Schwabe  & 
Scharmann  (1983)  and  Kamotani,  Ostrach  &  Vargas  (1983)  have  shown,  steady  convection  is  stable 


only  below  certain  values  of  the  Marangoni  number.  Above  the  critical  values,  the  flow  is  typically  an 
oscillatory  one,  with  commensurate  changes  in  the  transport  properties  and  important  implications  in 
material  processing  applications.  The  only  stability  analysis  of  this  class  of  flows  is  that  of  Smith  & 
Davis  (I983a,b)  who  restricted  their  attention  to  plane-parallel  flow  profiles  appropriate  to 
conduction-dominated  situations.  It  is  not  known  how  relevant  the  stability  limits  calculated  by  these 
authors  are  to  strongly  convective  situations. 

In  this  paper  we  compute  steady  thermocapillary  flows  in  a  square  cavity  (A  *  1)  by  a  finite- 
difference  procedure.  C  r  objective  is  to  obtain  accurate  numerical  solutions  to  the  stated  problem  at 
as  a  high  Reynolds  number  as  possible,  in  order  to  characterize  the  nature  of  strongly  convective 
flows  of  this  general  class.  Only  the  case  Ca  -»  0  is  considered,  so  that  the  free-surface  is  assumed 
flat  at  a  leading  order.  Surface  deflections  are  computed  by  domain  perturbation.  Boundary  layer 
formation  at  different  values  of  Pr  is  observed  at  large  values  of  Ma  and  Re.  We  use  the  numerical 
results  to  infer  the  relevant  scalings  of  a  consistent  boundary  layer  picture  of  the  flow,  valid 
asymptotically  as  Re  -»  oo.  We  do  not  attempt  to  present  a  complete  boundary  layer  theory  since 
solving  the  boundary  layer  problem  seems  to  be  as  difficult  as  solving  the  full  equations  of  motion. 


2.  Mathematical  model 


The  physical  model  consists  of  a  rectangular  cavity  of  average  height  d  and  width  w  containing  an 
incompressible,  Newtonian  liquid.  The  top  horizontal  boundary  is  a  free  surface  open  to  a  passive 
gas.  The  vertical  rigid,  isothermal  vails  are  differentially  heated  and  are  kept  at  temperatures  +  TQ/2 
relative  to  an  arbitrary  reference  temperature.  The  bottom  boundary  is  rigid  and  adiabatic.  In  the 
absence  of  gravity,  the  nondimensional  equations  for  the  liquid  motion  are  (Sen  &  Davis  1982) 


V-Y  -  0  . 

(2.1) 

Rev  -  (YY)  *  ~VP  +  VZY  . 

(2.2) 

Hav-(YT)  »  v2T  . 

(2.3) 

Here  length,  velocity,  temperature  and  pressure  are  dimensionless  with  respect  to  d,  yT Q//i,  TQ  and 
yTQ/d,  respectively,  where  /i  is  the  viscosity  and  surface  tension  is  assumed  to  decrease  with 
temperature  increase  at  a  constant  rate  y.  The  Reynolds  number  and  the  Marangoni  numbers  are 
defined  in  the  usual  way  by 

Re  =  yT0d/(/ir)  ,  (2.4) 

Ma  =  Re  Pr  .  (2.5) 

where  v  is  the  kinematic  viscosity  and  the  Prandtl  number  Pr  =  v  /  k,  and  k  is  the  thermal  diffusivity. 

The  motion  is  referred  to  a  Cartesian  coordinate  system  with  the  origin  at  the  middle  of  the  bottom 
boundary  with  the  y  axis  parallel  to  the  side  walls..  The  boundary  conditions  on  the  fixed  surfaces  are: 

Y(±l/(2A).y)  -  0  ;  K(x,0)  -  0  .  (2.6a,b) 

T(*1/(2A) ,y)  -  Tl/2  ;  Ty(x,0)  «  0  .  (2.7a,b) 

where  subscripts  denote  partial  differentiation  and  the  aspect  ratio  A  »  d  /  w  which  we  set  to  unity  in 
the  remainder  of  this  paper.  The  x  and  y  components  of  Y  will  be  denoted  as  usual  by  u  and  v, 
respectively. 

In  addition  to  Re,  Ma  and  Pr,  there  is  an  additional  dimensionless  parameter  which  is  a  measure  of 


the  free  surface  deformation.  This  is  the  capillary  number  Ca  and  is  given  by 


Ca  ■  yT0/<r0  .  (2.8) 

where  aQ  is  an  average  value  for  the  surface  tension.  We  see  that  Ca  ~  Ao/oQ,  and  in  experiments  is 
generally  a  small  quantity;  Kamotani  et.  al.  (1983).  We  thus  consider  the  case  Ca  — *  0,  i.e.  only  small 
variations  in  a  and  in  surface  deformation  are  allowed.  If  in  addition  to  small  Ca,  we  assume  the 
contact  angle  is  it/ 2,  the  surface  is  initially  located  at  a  height  y  =  1 .  Thus  to  leading  order  in  Ca,  we 
assume  a  flat  free  surface  with  the  boundary  conditions 

uy(x,l)  =  -Tx(x,l)  ;  Ty(x.l)  -  0  .  (2.9a.b) 

Corrections  to  the  surface  are  computed  as  follows.  If  we  denote  the  departure  of  the  free  surface 
from  y  » 1  by  h(x),  then  in  the  limit  Ca  — *  0  we  find 

hxx  "  ~Ca  P  f  h ( ±1/2 )  *  0  .  (2 . 10a, b) 

Equation  (2.10b)  assumes  a  fixed  contact  line.  Other  boundary  conditions  are  possible;  Sen  &  Davis 
(1982).  Because  the  pressure  field  p(x,y)  is  determined  only  to  an  additive  constant,  an  additional 
constraint  on  h(x)  from  global  continuity  is 

_1/2/1/2  h  dx  -  0  .  (2.11) 

Thus  h(x)  can  be  determined  to  0(Ca2)  by  solving  (2.10-2.11)  once  p  has  been  found.  We  note  for 
future  reference  that  the  surface  vorticity  «.  «  -ufx,l)  is  simply  the  surface  temperature  gradient. 


3.  Numerical  procedure 

Solutions  to  the  system  of  equations  (2.1),  (2.2),  (2.3),  (2.6),  (2.7)  and  (2.9)  must  be  constructed  by 
some  numerical  method.  At  a  first  glance  it  would  seem  that  a  spectral  method  with  an  infinite  order 
convergence  (Gottlieb  &  Orszag  1977)  may  be  suited  to  the  problem.  However,  the  hot  and  cold 
corners  on  y  =  1  are  singular.  The  vorticity 

«  =  vx  -  uy  .  (3.1) 

is  discontinuous  at  these  comers,  and  will  assume  different  values  as  a  corner  is  approached  on 
different  paths.  Indeed,  following  the  analysis  of  Moffatt  (1964)  it  can  be  shown  that,  in  the  vicinity  of 
the  hot  corner,  the  vorticity  is  given  by 

w  *  -Tx(-1/2,1)  (1-4 e/v),  (3.2) 

where  the  angle  6  is  0  on  x  »  -1/2  and  v/2  on  y  =  1.  It  should  be  noted  that  the  flow  field 
corresponding  to  (3.2)  is  valid  only  very  close  to  the  comer  so  that  the  flow  is  locally  a  Stokes  flow.  A 
similar  expression  holds  for  u  near  the  cold  comer.  It  is  thus  expected  that  <j  will  vanish  on  the 
diagonals  as  the  comers  are  approached. 

Since  this  discontinuity  will  limit  the  convergence  of  any  spectral  method  we  opt  for  a  finite 
difference  approach.  The  method  and  the  iteration  procedure  we  use  to  solve  the  finite  difference 
equations  are  described  in  detail  by  Patankar  (1981).  Briefly,  the  computational  region  is  divided  into 
rectangular  control  volumes  with  the  grid  points  located  at  the  geometric  centers  of  these  cells. 
Additional  boundary  grid  points  are  included  where  the  boundary  conditions  (2.6, 7, 9)  are  imposed. 
The  finite  difference  equations  are  obtained  by  integrating  the  governing  equations  (2.1, 2, 3)  over  the 
control  volumes  with  assumed  local  linear  variations  in  any  of  the  primitive  variables.  The  convection 
and  diffusion  fluxes  are  approximated  by  a  power-law  scheme.  The  staggered  location  for  the 
velocity  components  is  adopted  and  the  velocity- pressure  coupling  treatment  follows  Patankar’s 
(1981)  revised  procedure.  A  line-by-line  iteration  to  solve  the  discretized  equations  is  used  with  one 
completed  iteration  comprising  five  double  sweeps  of  the  field.  Under- relaxation  in  solving  for  u,  v 
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and  T  was  required;  a  relaxation  factor  of  0.85  was  used  throughout. 

After  a  number  of  solutions  were  computed  on  uniform  grids  with  65  x  65  and  80  x  80  mesh  points, 
we  used  a  nonuniform  grid  with  62  points  in  the  x  direction  and  54  points  in  the  y  direction,  with  the 
highest  concentration  of  the  mesh  points  near  the  cold  comer  where,  as  we  shall  see,  the  sharpest 
gradients  occur.  The  smallest  control  volume  at  the  cold  comer  is  a  square  with  side  length  0.005. 
The  rectangular  cell  at  the  hot  corner  has  an  x  extent  of  0.01 .  The  mesh  spacing  was  gradually 
increased  away  from  the  cold  and  hot  corners  in  both  the  x  and  y  directions.  The  largest  mesh 
spacing  near  the  middle  of  the  bottom  boundary  is  less  than  0.05.  All  computations  were  performed 
on  a  VAX  11/780.  A  converged  solution  requires  from  250  to  450  iterations,  with  each  iteration 
completed  in  less  than  25  CPU  seconds.  Convergence  was  assumed  when  the  largest  variation  in  any 
of  u,  v,  p  and  T  was  less  than  some  convergence  tolerance  which  we  set  to  10'5.  For  extreme  values 
of  Pr,  we  also  performed  some  computations  on  a  nonuniform  mesh  with  the  same  smallest  elements 
but  with  70  x  60  mesh  points. 

In  table  I  we  list  some  representative  computed  results  with  different  grids  from  which  we  claim 
accuracy  within  3%  for  results  with  the  62  x  54  grid.  It  should  be  noted  that  the  largest  "error’'  is  in 
the  magnitude  of  of  circulation  •4fmax  (here  <p  is  the  streamfunction).  This  is  especially  true  at  large 
Re  and  when  the  point  at  which  occurs  is  located  far  from  the  top  comers,  i.e.  where  mesh 
spacing  is  relatively  large.  Because  we  use  the  conservative  (divergence)  form  of  the  governing 
equations,  an  exact  solution  of  the  difference  equations  should  result  in  equal  values  for  the  Nusselt 
numbers  computed  on  the  hot  and  cold  boundaries  Nu  and  Nu+,  respectively.  Decreasing  the 
convergence  tolerance  results  in  yet  closer  values  for  Nu  and  Nu  +  than  those  listed  in  table  I. 


4.  Numerical  results;  Pr  s  1 

Our  primary  interest  is  in  the  character  of  the  motion  at  large  values  of  Ma  and  Re.  Equation  (2.3) 
shows  that  Ma  measures  the  strength  of  temperature  (thermal  energy)  convection  to  diffusion.  Thus 
large  values  of  Ma  will  lead  to  the  formation  of  thermal  boundary  layers.  Equations  (2.1)  and  (2.2)  can 
be  reduced  to  the  vorticity  transport  equation 

Rey.(Vw)  *  v2«  .  (4.1) 

and  hence  large  Re  implies  the  formation  of  regions  of  concentrated  vorticity,  in  particular  near  the 
free  surface.  We  first  consider  the  case  Pr  >  1,  i.e.  simultaneous  formation  of  both  temperature  and 
vorticity  layers  with  Re  =  Ma. 

With  Pr  *  1,  we  compute  the  thermocapillary  motion  for  102<  Re<  104.  As  expected,  the  flow 
consists  of  a  major  closed  circulation,  accompanied  by  Moffatt  comer  eddies  at  the  two  lower 
corners.  The  temperature  is  distorted  from  a  linear  conductive  field  by  this  circulation,  especially  at 
high  Ma.  The  surface  temperature  distribution  is  of  key  interest,  since  it  is  the  gradient  of  this  quantity 
which,  through  the  thermocapillary  stress,  drives  the  motion.  As  an  example  of  these  features,  we 
show  in  Figure  la-c  the  streamlines,  isotherms  and  iso- vorticity  contours  for  Re  =  Ma  =  103.  As  can 
be  seen  from  Figure  1c,  there  is  a  large  concentration  of  vorticity  near  the  stagnation  point  near  the 
cold  (right)  boundary.  It  is  further  seen  that  the  stagnation  point  near  the  hot  (left)  boundary  exhibits 
little,  if,  any,  boundary  layer  behavior.  The  dividing  u  =  0  line  is  oriented  at  46°,  as  required  by  (3.2). 

In  Figure  2  we  show  the  results  of  our  parametric  study  at  Pr  =  1.  We  plot  surface  temperature 
T(x,l),  velocity  u(x,1),  and  the  surface  deflection  h(x)  corresponding  to  five  different  values  of  Re  (or 
Ma).  For  Re  less  than  about  1000,  it  is  seen  from  Figure  2a  that  the  surface  temperature  gradient  at 
the  hot  comer  Tx(-1/2,1)  is  greater  than  the  conduction  value  of  -1,  i.e  a  decreased  surface  vorticity 
at  the  hot  comer.  At  these  relatively  low  values  of  Re,  Figure  2b  indicates  that  the  surface  velocity 
increases  monotonically  from  zero  and  drops  back  smoothly  to  zero  at  the  cold  corner.  The 
associated  surface  deflection  plot  in  Figure  2c  shows  that  there  are  two  peaks  in  h(x)  which  decrease 
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in  magnitude  and  move  towards  the  corners  with  increasing  Re.  With  Re  (or  Ma)  increasing  above 
1000  however,  Tx(-1/2,1)  is  less  than  -1  and  decreasing,  u(x,1)  exhibits  two  peaks  and  h(x)  showing  a 
more  or  less  flat  depression  on  an  increasing  length  of  the  free  surface  near  the  hot  (left)  boundary.  It 
is  also  observed  that  the  surface  elevation  near  the  cold  corner  is  increasing  with  Re.  From  these 
features  it  may  be  concluded  that  a  boundary  layer  regime  is  reached  for  Re  (  or  Ma)  greater  than 
about  1000. 


The  streamlines,  vorticity  and  temperature  fields  corresponding  to  Re  =  Ma  =  104  are  shown  in 
Figure  3.  From  this  Figure  it  may  be  concluded  that  the  limiting  flow  field  at  large  values  of  Re  (or  Ma) 
will  consist  of  an  irrotational,  isothermal  core  with  boundary  layers  forming  on  the  four  boundaries. 
The  corner  eddies  predicted  by  Moffatt  (1964)  are  present  in  Figure  3a,  while  the  vorticity  is 
discontinuous  near  the  corners  in  accord  with  equation  (3.2).  The  sharpest  gradients  occur  near  the 
cold  comer  as  evident  in  Figure  3c. 


Of  interest  is  the  total  heat  transport,  or  the  Nusselt  number,  defined  by 
Nu*  *  o/1  -Tx(*l/2.y)  dy  , 


(4.2) 


As  noted  above,  Nu+  =  Nu  to  within  1%.  In  Figure  4  we  show  the  variation  of  Tx  along  the  hot  and 
cold  walls  from  which  it  follows  that  while  all  of  the  hot  boundary  participates  in  the  heat  transfer 
process,  only  a  "small"  cold  corner  region  contributes  to  most  of  the  transfer.  This  picture  will  be 
validated  in  section  5. 


The  vigor  of  the  convective  motion  is  indicated  by  the  magnitudes  of  the  heat  transfer,  flow 
circulation  ('/'max)  and  maximum  vorticity  (w^).  Figure  5  is  a  logarithmic  plot  of  the  variation  of 
Nu  *  (Nu  +Nu+)/2,  •<>>wax  and  4  with  Re,  where  u  *  Tx(1/2,1)  for  all  the  cases  with  Pr  =  1 


that  we  computed.  From  Figure  5  it  may  be  concluded  that 


0(Re2'3)  , 


Hu  ~  0(Re1/3)  . 


(4.3) 

(4.4) 
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5.  Structure  of  boundary  layers;  Pr  =  1 
Below  we  develop  a  self-consistent  picture  of  the  boundary  layer  structure  which  exists 
asymptotically  as  Re  -+  oo.  We  show  that,  beginning  with  the  key  assumption  that  surface  vorticity 
becomes  independent  of  Re  as  Re  -*  oo,  all  of  the  available  numerical  evidence  is  in  agreement  with 
the  inferred  scalings.  The  first  region  of  interest  is  the  boundary  layer  near  the  free  surface. 

Free  surface  layer 


Here  x  ~  0(1).  The  assumption  u  —  0(1),  which  follows  from  the  boundary  condition  (2.9a)  and 
Figure  2c,  together  with  the  required  balance  of  diffusion  and  convection  lead  to  the  free  surface 
transformation: 


y  =  1  -  Re'1/3y 


(S.la.b) 


X  =  X 

u  =  Re'1/3u(x,y)  ;  v  =  Re‘2/3v( x ,y )  ;  * p  a  Re'2/3J(x,y)  ,  (5.1c,d,e) 

p  =  Re1/3p(x,y) 

T  *  T(I.y)  . 


(5. If) 

(5-lfl) 


Equation  (5.1c)  implies  a  surface  velocity  decreasing  with  Re  in  qualitative  agreement  with  Figure  2b. 
The  boundary  layer  equations,  to  leading  order  become 


u-  +  v-  *  0 

x  y 


uu,  -  vuy  *  -Ps  ♦  5y- 
0  =  pv  . 


uT.  -  It 


7 


Pr  ‘TJy  • 


(5.2a) 
(5.2b) 
(5.2C) 
(5. 2d) 


These  are  the  usual  boundary  layer  equations.  The  known  boundary  conditions  are 


y  3  0  :  v  =  ?-  =0  ;  -u-  »  T-  . 


(5.3a,b,c) 


Matching  conditions  for  y  -»  oo  with  a  core  solution  and  initial  conditions  for  x  -»  -1/2  must  be 
provided. 


Core  region 


Formally,  the  limiting  equations  in  the  core  are 

(Y.y)w  =  0  ,  (5.4a) 

(Y.V)T  3  0  .  (5.4b) 

Thus  ,  in  principle,  u  and  T  are  constants  along  streamlines.  Figures  2b  and  2c  indicate  that  the 
limiting  core  flow  is  isothermal  (at  some  Tc)  and  irrotational.  The  pressure  field  is  then  found  from  the 
Bernoulli  equation 

p  +  Re  Vz/2  *  constant  .  (5.5) 

Matching  with  the  free  surface  layer  indicates  that 

p  ~  0(Re1/3)  ;  V  ~  0(Re‘1/3)  .  (5.6a,b) 

so  that  the  core  streamfunction 

V-  ~  0(Re‘1/3)  .  (6.7) 

This  is  consistent  with  equation  (4.6)  and  Figure  4. 

Cold  comer  region 

This  stagnation  point  flow  region  is  where  most  of  the  heat  transfer  on  the  cold  boundary  occurs 
(Figures  2a,  3c,  4b)  and  where  the  largest  vorticity  is  generated.  The  horizontal  free  surface  flow 
turns  in  this  region  as  u(x,1)  drops  sharply  to  zero  (Figure  2b).  The  y  extent  and  the  amount  of  flow  in 
this  region  are  determined  by  those  in  the  free  surface  region.  In  addition,  the  required  convection, 
diffusion  balance  leads  to  the  cold  comer  scalings 


y  «  1  -  Re'1/3  y 

;  x  ■  1/2  -  Re*z/3  x  , 

(5.8a,b) 

u  *  Re*1/3  u(x,y) 

;  v  ■  v(x,y)  ;  ^  «  Re'z/3  $(x,y) 

(5.8c,d,e) 

u  ■  Re2/3  w(x.Jf) 

• 

(5.8f) 

p  -  Re1'3  p(S,y) 

• 

(5.8g) 

T  »  T(x.y)  . 

(5.8h) 
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The  leading  order  equations,  boundary  and  matching  conditions  become 


u7  +  v-  -  0 

-UU;  +  vu.  =  p;  +  uJx-  . 

-Cv.  -  VV-  *  V~  . 

-uf.  -  vf-  =  Pr'1?;;  . 

x  =  0  :  u  =  v  »  0  ;  T  =  -1/2  . 

x  — »  oo  :  v  -»  0  ,  u(x-*oo,y)  —  u(x-*l/2,y)  . 


(5.9a) 

(5.9b) 

(5.9c) 

(5.9d) 

(5.9e.f.g) 

(5.9h.i) 


These  quantities  can  not  satisfy  the  initial  conditions  on  y  ■  0  since  u-  ~  0(1)  while  T-  —  0(Re2/3). 

y  * 

Thus,  there  will  be  a  corner  subregion  where 


1/2  -  x  ~ 

0(Re'2/3)  . 

(5.10a) 

1  -  y  ~ 

0( Re'1)  . 

(5.10b) 

*  ~ 

0(Re*4/3)  . 

(5.10c) 

T  ~ 

0(1) 

( 5 . lOd) 

It  should  be  noted  that  (5.10c)  is  a  consequence  of  higher  order  free  surface  layer,  while  (5.10b)  is 
necessary  to  satisfy  the  vorticity  balance  on  the  free  surface. 


Equations  (5.8)  is  in  complete  agreement  with  Figures  3b,  4,  indicating  that  most  of  the  heat  transfer 
on  the  cold  wall  is  within  the  corner  region.  In  addition,  it  follows  from  equation  (4.3)  that 

Nu  *  „/>  -Tx( 1/2 ,y)  dy  . 

~  Re 1/3  0/°°  T;(0.y)  dy  .  (5.11) 

which,  again,  is  in  agreement  with  Figure  5. 


The  existence  of  these  regions  with  various  asymptotic  scales  (see  equations  5.1, 5.8  and  5.10),  is 
very  demanding  on  the  computational  mesh.  This  is  particularly  so  at  extreme  values  of  Pr. 
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6.  Numerical  results;  Pr  i  1 

Equations  (2.2)  and  (4.1)  show  that  convection  of  vorticity  is  stronger  (weaker)  than  convection  of 
energy  according  to  Pr<1  (Pr>1).  The  influence  of  Pr  on  the  motion  is  ascertained  from  a  sequence  of 
computations  with  Pr  =  0.05,  0.1,  1,  10  and  50.  The  largest  value  of  Re  or  Ma  attainable  with 
reasonable  accuracy  was  50,000. 

Figures  7a,b  are  plots  of  the  surface  temperature  corresponding  to  Pr  *  0.1  and  50  at  various 
values  of  Re.  It  is  seen  that,  with  increasing  Re,  the  surface  vorticity  at  the  hot  corner  first  decreases, 
and  then  begins  to  increase  monotonically  with  further  increase  in  Re.  This  is  similar  to  the  case  Pr  » 

1  in  Figure  2a.  The  surface  velocity  corresponding  to  the  parameter  values  of  Figures  7a, b  are  shown 
in  Figures  8a, b.  It  is  observed  that  a  two-peak  structure  is  eventually  approached  at  sufficiently  large 
Re  (or  Ma).  It  is  important  to  note;  Figure  8b,  that  at  appropriately  low  values  of  Re  (depending  on  Pr) 
that  for  Pr  >  1  the  motion  is  faster  near  the  hot  corner  than  almost  every  where  on  the  free  surface 
except  in  the  readily  developed  thermal  (and  vorticity  layer  at  the  cold  comer).  The  opposite  is  true 
for  Pr  <  1  and  low  Re;  Figure  8a,  where  convection  is  expected  to  dominate  near  the  cold  comer. 

The  influence  of  Pr  on  the  pattern  of  circulation  is  evident  from  Figures  9a, b  where  show  the 
streamlines  of  the  thermocapillary  motion  corresponding  to  Pr,  Re  and  Ma  values  of  0.05;  1000;  50 
and  50;  200;  10,000.  At  sufficiently  low  Re  it  is  seen  that  the  point  where  occurs  is  close  to  the 
cold  comer  when  Pr  <  1,  while  for  Pr  >  1  this  point  occurs  near  the  hot  comer.  Thus,  again,  it  is 
concluded  that  convective  effects  are  more  important  at  the  hot  (cold)  corner  according  to  Pr>  1  (<  1) 
at  appropriate  low  values  of  Re.  With  increasing  Re,  however,  the  pattern  of  motion  is  similar  to  that 
for  Pr  «  1  in  Figure  3a.  This  behavior  with  Pr  was  also  found  by  Fu  &  Ostrach  (1983)  in  their 
axisymmetric  half-zone  model. 

Plots  of  •'Pmax,  Nu  and  versus  Re,  similar  to  those  in  Figures  5a,b,c  for  Pr  *»  1,  show  that  the 
asymptotic  estimates  given  in  equations  (4.4),  (4.5)  and  (4.6)  are  applicable  with  Pr  ■  1.  This  is  also 
indicative  that  the  character  of  the  motion  becomes  independent  of  Pr  at  sufficiently  large  Re. 


•  -  •  «  *  »  ’  •  •  *  s  .»  v  v  v  '_**  *  *  V 


The  most  remarkable  influence  of  Pr  on  thermocapillary  convection  is  found  on  the  shape  of  the 
free  surface.  Figures  lOa.b  show  h(x)  for  the  parameter  values  of  Figures  7a,b  and  8a,b.  It  is  seen 
that  for  Pr  >  1,  the  region  of  strong  motion  near  the  hot  corner  is  accompanied  by  sufficient  low 
pressure  that  the  largest  depression  exceeds  the  largest  elevation.  At  Pr  <  1  however,  there  is  a  build 
up  of  pressure  sufficient  to  produce  a  secondary  elevation  near  the  hot  corner. 

It  is  of  interest  to  note  that,  for  all  the  cases  considered,  |h(x)|  is  small,  (Figures  2c,  10),  indicating 
that  the  range  of  capillary  numbers  for  which  the  surface  deflection  is  accurately  given  by 
perturbation  in  Ca  may  be  rather  large. 


■  -  W  <■  v  ‘.'f-.*  v  v 
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7.  Concluding  remarks 

We  have  presented  the  results  of  a  reasonably  complete  study  of  thermocapillary  convection  in  a 
square  cavity.  Our  accurate  computational  procedure  allows  us  to  consider  situations  with  large  Re 
and  Ma,  and  thus  observe  the  formation  of  boundary  layers,  in  particular  at  the  cold  stagnation  point 
region.  This  boundary  layer  structure  is  shown  to  be  consistent  with  an  asymptotic  theory  valid  as 
Re— ♦  oo. 

We  encountered  no  difficulty  in  computing  two-dimensional  steady  states  by  time-like  iterations.  It 
is  also  felt,  from  numerical  experiments,  that  such  steady  motions  continue  to  exist  at  yet  higher 
values  of  Re  and  Ma.  Since,  from  previous  experimental  work  and  from  the  analysis  of  Smith  &  Davis 
(1983a),  it  is  expected  that  oscillatory  motions  will  prevail  at  some  “supercritical"  Re  and  Ma,  it  is 
conjectured  that  such  an  unsteady  motion  is  three-dimensional.  It  seems  that  there  is  no  convenient 
method  to  study  the  instability  of  the  solutions  we  computed  on  the  finite  difference  grid  to  three- 
dimensional  disturbances.  Thus,  further  progress  in  the  study  of  this  phenomena  must  be  both 
three-dimensional  and  time  dependent. 
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Influence  of  finite-difference  grid  on  the  solution 


Pr 

Re 

M  x 

N* 

Nu_ 

Nu, 

-\p  .  x  102 

*u 

■4 

max 

0.05 

5  x  104 

62  x 

54 

2.120 

2.131 

0.174 

14.37 

70  x 

60 

2.161 

2.151 

0.180 

14.44 

0.1 

5  x  104 

62  x 

54 

2.937 

2.957 

0.168 

29.04 

70  x 

60 

2.980 

2.981 

0.175 

28.95 

1 

5  x  103 

65  x 

65 

3.459 

3.448 

0.370 

28.72 

80  x 

80 

3.454 

3.447 

0.373 

31.20 

62  x 

54 

3.420 

3.412 

0.366 

38.37 

50 

5  x  10z 

62  x 

54 

4.895 

4.898 

0.155 

139.4 

70  x 

60 

4.894 

4.896 

0.155 

139.7 

•  M  and  N  are  the  number  of  grid  points  in  the  x  and  y  directions  respectively.  The  65  x  65  and  80  x 
grids  are  uniform  while  the  62  x  54  and  the  70  x  60  grids  are  graded. 
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Figure  captions 

•  Figure  1.  Thermocapillary  flow  at  Pr  =  1  and  Re  =  1000.  a)  Streamlines  at  equal 

increments  of  circulation.  The  motion  is  clockwise  with  4  =  0.0048.  b)  Isotherms  at 

equal  increments,  c)  Iso-vorticity  contours.  u  =  -1.0  and -11.79  at  the  top  hot  (left)  and 
cold  (right)  corners,  respectively.  The  largest  positive  u  is  7.3. 

•  Figure  2.  a)  Surface  temperature  corresponding  to  Pr  =  1  and  Re  a  100,  500,  1,000, 
5,000  and  10,000.  b)  Associated  surface  velocity,  and  c)  surface  deflection. 

•  Figure  3.  Same  as  Figure  1  but  with  Re  =  10,000.  Here  4max  =  0.003;  u  =  -2.4  and 
•60.2  at  the  top  hot  and  cold  corners,  respectively;  and  the  largest  value  of  u  is  10.5. 

•  Figure  4.  Variation  with  y  of  Tx  on  the  a)  hot  boundary,  and  b)  cold  boundary  at  Pr  a  i 
and  Re  a  1,000, 5,000  and  10,000. 

•  Figure  5.  Variation  with  log  Re  of  log  (-</'),  log  (-*»  >  and  log  Nu.  Also  shown  are 

max  max 

straight  lines  with  slopes  -1/3, 2/3  and  1  /3. 

•  Figure  6.  Sketch  of  the  important  surface  layer  regions.  The  following  are  the  scalings  in 
regions  I,  II,  III  and  the  core. 
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1 
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Re2'3 
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Re'1/3 

Re'1'3 

Re*1'3 

1 

1 

•  Figure  7.  Surface  temperature  T(x,1)  at  a)  Pr  =  0.1  and  Re  =  103,  104,  2  x  104  and 
5  x  104.  b)  Pr  a  50  and  Re  a  20, 100,  200  and  500. 

•  Figures  8a, b.  Associated  surface  velocity  u(x,l)  at  the  parameter  values  of  Figures  7a, b, 
respectively. 

•  Figure  9.  Streamlines  corresponding  to  a)  Pr  =  0.05  and  Re  =  1000;  4max  *  0.0088,  and 
b)  Pr  -  50  and  Re  -  200;  4^  »  0.0019. 

•  Figure  10a, b.  Surface  deflection  at  the  parameter  values  of  Figures  7a, b,  respectively. 
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